Swarm behavior of self-propelled rods and swimming flagella 
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Systems of self-propelled particles are known for their tendency to aggregate and to display swarm behav- 
ior. We investigate two model systems, self-propelled rods interacting via volume exclusion, and sinusoidally- 
beating flagella embedded in a fluid with hydrodynamic interactions. In the flagella system, beating frequencies 
are Gaussian distributed with a non-zero average. These systems are studied by Brownian-dynamics simulations 
and by mesoscale hydrodynamics simulations, respectively. The clustering behavior is analyzed as the particle 
density and the environmental or internal noise are varied. By distinguishing three types of cluster-size prob- 
ability density functions, we obtain a phase diagram of different swarm behaviors. The properties of clusters, 
such as their configuration, lifetime and average size are analyzed. We find that the swarm behavior of the two 
systems, characterized by several effective power laws, is very similar. However, a more careful analysis re- 
veals several differences. Clusters of self-propelled rods form due to partially blocked forward motion, and are 
therefore typically wedge-shaped. At higher rod density and low noise, a giant mobile cluster appears, in which 
most rods are mostly oriented towards the center. In contrast, flagella become hydrodynamic ally synchronized 
and attract each other; their clusters are therefore more elongated. Furthermore, the lifetime of flagella clusters 
decays more quickly with cluster size than of rod clusters. 



I. INTRODUCTION 

Systems of self-propelled particles (SPP), which exhibit 
an interaction mechanism that favors velocity alignment of 
neighboring particles, often display collective behaviors like 
swarming and clustering. There are many examples for this 
swarming behavior, ranging from systems of microscopic par- 
ticles (sperm, bacteria, nano-rods) to systems of macroscopic 
objects (birds, fish). 

Since the pioneering simulation work of Vicsek et al. (TJ, 
SPP systems have attracted a lot of interest at the theoreti- 
cal l|2T[8l and computational ll9i fT5l level. Typically, in sim- 
ulation models of swarm behavior, point-like agents move 
with an imposed non-zero velocity and tend to align their di- 
rection of motion with others in a prescribed neighborhood 
|[T][T0l[TTl[T4). Although the alignment mechanism may differ 
from one model to the other, the basic properties of swarm be- 
havior are quite universal 1 16 1. Upon variation of parameters 
such as particle density, particle velocity, or environmental 
noise, the system can undergo a transition from a disordered 
state, where the average total velocity or orientation vanishes, 
to a nematically ordered state. Near the transition point, the 
cluster-size probability density function is characterized by a 
power-law decay ifTTl [TBI . For intermediate densities, phase 
separation into regions of different density and band formation 
has been found lfT31l . 

Self-propelled motion is common in biological systems at 
micro- or mesoscopic length scales, such as suspensions of 
bacteria, like E. coli [17] and Bacillus subtilis lfT8l - E0l . or tis- 
sue cells (keratocytes) [9|, whose sizes are all on the microm- 
eters scale. A special class of biological systems are rod-like 
self-propelled particles (rSPP), for example myxobacteria (ap- 
proximately 10;um long) ll2Tll22ll . When starved, myxobacte- 
ria are elongated to an average aspect ratio of approximately 
1:7, glide on a substrate along their long axis and undergo 
a process of alignment, rippling, streaming and aggregation 
that culminates in a three-dimensional fruiting body. A model, 



which takes into account the exchange of a morphogen during 
cell-cell contact and a preferred cell motion in the direction 
of largest morphogen concentration, has been designed to de- 
scribe the streaming and two-stage aggregation of myxobac- 
teria ED . 

Sperm (with a length of about 50^m) l24l l25ll and nema- 
todes [26 1 (about 1mm long) employ a sinusoidal undulation 
of their slender bodies to push the fluid backwards and to 
propel themselves forward. Large train-like clusters of wood 
mouse sperm l27l l28l are believed to result in greater thrust 
forces to move more efficiently through a highly viscous envi- 
ronment. The wood mouse sperm has a hook-like structure at 
its head, by which it can be hitched to the mid-part or the tail 
of a neighboring cell for robust cooperation. However, nema- 
todes which do not have hook structures, also display a pro- 
nounced tendency to adhere to each other in a film of water, 
to form assemblies consisting of many organisms, and to ex- 
hibit a striking co-ordinated movement [26 1 . Also, sea urchin 
sperm organize into a hexagonal pattern of rotating vortices at 
surfaces 

A nice physical realization of self-propelled rods (SPR) are 
bimetallic nano-rods consisting of long Pt and Au segments 
OUll . The rods, about 300«m in diameter and 2fim long, move 
autonomously in an aqueous hydrogen peroxide solutions by 
catalyzing the formation of oxygen at the Pt end. They move 
predominantly in the direction of the Pt end, with a velocity 
depending on the concentration of hydrogen peroxide. When 
a gradient of the hydrogen peroxide concentration is imposed, 
the rods exhibit directed motion towards regions of higher 
concentrations through active diffusion f3"Tl . 

A related system is a fluidized monolayer of macroscopic 
rods in the nematic liquid crystalline phase [32 1 . The rods 
confined between two hard walls are energized by an exter- 
nal vertical vibration, and gain kinetic energy through fre- 
quent collisions with the floor and the ceiling of the container. 
Long-lived giant number fluctuations are found, which shows 
that simple contact can give rise to flocking, coherent swirling 
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motion and large-scale inhomogeneities [33]. However, in this 
experiment, the rods do not have a preferred direction of mo- 
tion. 

All of these examples of self-propelled particles employ 
different propulsion mechanisms and have different interac- 
tions. However, their swarm behavior, such as flocking, 
streaming and clustering, is surprisingly similar. The com- 
mon characteristic of these systems is their rod-like structures 
and their quasi-two-dimensional active motion. Myxobacteria 
glide on surfaces BTTl . while sperm and nematodes gather at 
substrates [26, 29 34 1. In suspensions of rod-like particles in 
thermal equilibrium, volume exclusion favors the alignment 
of rods. At high densities, it stabilizes a nematic state charac- 
terized by long-range orientational order [35|. 

While constant- velocity polar point particles interacting lo- 
cally by nematic alignment in the presence of noise have 
been studied intensively in recent years |[T HT5ll . much less is 
known theoretically about the behavior of elongated particles 
with volume exclusion, or about the collective behavior of 
swimmers with hydrodynamic interactions. Previous simula- 
tion studies of self-propelled rods (SPR) in two dimensions 
show that self-propelled motion enhances the tendency for 
nematic ordering P6l . as well as aggregation and clustering 
||371 . Also, rods have an increased probability to be located 
near surfaces (depending on their velocity, length and ther- 
mal noise) 1 38 1 and form hedgehog-like clusters at surfaces 
11391 . In Ref. IT371 . two regimes of clustering have be distin- 
guished by their unimodal or bimodal weighted cluster-size 
distribution functions; however, the system contained a rela- 
tively small number of particles compared to those employed 
in simulation studies of swarming of SPPs. Continuum equa- 
tions for the description of SPR systems have been derived re- 
cently within a mean-field approximation ||6l |7) . This theory 
predicts that hard-core interactions are insufficient to generate 
a macroscopically polarized state, because they cannot distin- 
guish the two ends of a rod, and makes interesting predictions 
for the fluctuations in the nematic and isotropic state (such 
as a crossover from diffusive to propagating density fluctua- 
tions). However, the mean-field approximation of volume ex- 
clusion has the limitation of omitting correlation effects, and 
thus works best for slowly varying density distributions. 

In addition, hydrodynamic interactions between rSPP have 
so far been largely neglected. These interactions depend on 
the type of self-propulsion, where "pullers" repel and "push- 
ers" attract each other I40ll4ll . Nematic suspensions of swim- 
ming rod-like pushers are found to be unstable at long wave- 
lengths as a result of hydrodynamic fluctuations B2l . For 
sperm and flagella, it has been shown theoretically that the 
hydrodynamic coupling synchronizes the phases of their si- 
nusoidal beating tails 11241 [43l l44ll . Also, the hydrodynamic 
interaction between these microswimmers implies attraction 
and cluster formation ll43l : similarly, it makes an essential 
contribution to the capturing of sperm near walls [45 1. How- 
ever, the relative importance of directed self-propulsion, par- 
ticle shape, volume exclusion, and hydrodynamic interactions 
to the emergence of swarm behavior remains unclear. 

In this paper, we employ a model of hard rods with strict 
volume exclusions and simulate large systems containing at 



least 1000 particles. We focus on rSPP systems at a density 
below the isotropic-nematic transition of Brownian rods. We 
employ a model consisting of rigid SPR performing an over- 
damped translational motion in two dimensions, and analyze 
the resulting cluster-size probability density distribution, clus- 
ter configurations and lifetimes. Three types of cluster-size 
probability density distribution functions allow to distinguish 
three different states, and to construct a phase diagram as a 
function of particle density and environmental noise. As a 
special case of rSPP with an explicit propulsion mechanism, 
we investigate a suspension of flagella, which move by sinu- 
soidal beating of their body in a two-dimensional fluid. The 
motion of the surrounding fluid is described by particle-based 
mesoscopic simulation method called multi-particle collision 
dynamics (MPC) ll46l l47l . This method has been shown to 
capture the full hydrodynamics and flow behavior of complex 
fluids over a wide range of Reynolds numbers very well [4-8 1 . 
By comparing the results for SPR and flagella, we elucidate 
the contribution of hydrodynamic interactions to the swarm 
behavior. 

This paper is organized as follows. Section [II] gives a brief 
description of our models and simulation methods. We an- 
alyze the collective behavior of SPR systems in Sec. 
Sec. 
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we study the swarm behavior of flagella, and com- 
pare the results obtained with both models. The influence of 
hydrodynamic interactions and the flagellar beat on the swarm 
behavior are discussed. We summarize our main conclusions 
in Sec.[V] 



II. MODELS AND SIMULATION TECHNIQUES 

A. Self-Propelled Rods 

We consider a system of N ro d rods of length L ro d in a two- 
dimensional simulation box of size L x x L y . Each rod is char- 
acterized by an orientation angle 6 ro d,i with respect to the x- 
axis, a center-of-mass position r ro dj, a center-of-mass velocity 
\ ro d.i and an angular velocity (o r od,i around its center of mass 
(see Fig. The rods move ballistically according to their 
velocities, 

r rodJ (t + At ro d) = r TO d,i{t) + Vrod,i(t) At rod , (1) 
6 ra d,i(t + Atrod) = OrodA*) + Urod,i(t)At ro d , (2) 

where At ro d is the simulation time step. The particle velocity 
can be decomposed into a parallel and a perpendicular com- 
ponent relative to the rod axis, v r od,i = v rod,i,\\ + v rod,i,±- 

We consider the rods to be embedded in an overdamped 
fluid medium where hydrodynamics can be approximated by 
an anisotropic friction on the rod-like particles. The motion is 
then determined by 
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FIG. 1: (Color online) (a) Model of a self-propelled rod, and the 
coordinates used in two dimensions. The rod is discretized into 
n ro< j = Lyodlh beads for the calculation the volume-exclusion inter- 
action, (b) Model of a flagellum in two dimensions. 
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where eg and e ± are the local parallel and perpendicular unit 
vectors of the rod orientation. F ro d,o is a constant propelling 
force applied along ey. The friction coefficients are given by 
y ± = 2y\\, 7n = L r „d and y r = y\\L ro d 2 /6. The random forces 
fll, £l and % r are white noises, which are are determined by 
their variances <T 2 od L rod , o- 2 rod L rod and v^L^J 12, respec- 
tively. Finally, F, y is the force generated by volume exclusion 
between rods i and j, and My is the torque generated by F, y 
on rod i in the reference system of center of mass of rod i. 

For the calculation of the interactions, each rod is dis- 
cretized into n w d - L ro dllb beads of diameter l b , as illustrated 
in Fig. [ij. The volume exclusion between rods is then mod- 
elled by a shifted and truncated Lennard- Jones potential 
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between beads belonging to different rods, where r is the dis- 
tance between two beads, l b is the bead diameter, and e is the 
strength of the potential. We use e as the energy scale in our 
SPR simulations. 

A single rod without noise then moves with a constant 
velocity vo = F ro dply\\. In the non-zero noise regime, 
the diffusion constant along the parallel direction is D\\ - 
cr 2 rgd L ro dAt ro d /2y 2 . The dimensionless Peclet number, which 
measures the ratio of self-propelled and diffusive motion, is 
thus 



Pe = 



cr rod Mrod 
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We use 1 /Pe oc cr 2 od to characterize the strength of the envi- 
ronmental noise. 

In SPR systems 10 [37), alignment is naturally introduced 
by the volume exclusion between the anisotropic particles; 
this also implies that the interaction neighborhood needs no 
further assumptions, but is directly related to the rod length. 
Hard-core interactions do not distinguish the two ends of an 



symmetrically elongated object. Thus, both parallel and anti- 
parallel velocity configurations are induced. In simulations 
of point-like SPPs, noise is implemented by adding a random 
component to the velocity orientation of each particle. In our 
model of SPR, random forces are applied on each rod, which 
results in fluctuations in both the magnitude and the orienta- 
tion of the velocity vectors. For a single rod, the orientation 
fluctuations lead to rotational diffusion, which implies a per- 
sistence length 
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of its trajectory. Note that the noise forces are not caused 
by thermal fluctuations, which would require a factor two 
between the variance of the random forces in parallel and 
perpendicular directions. In most biological and synthetic 
rSPP systems, thermal fluctuations are indeed negligible due 
to large size of the particles. In these systems, the environ- 
mental noise arises, for example, from density fluctuations of 
signalling molecules for chemotactic swimmers, or from fluc- 
tuations of the motor activity. 

We use rods of length L ro d — 1 \l b and undisturbed velocity 
Vo = l.2le/(yL rod ). Effects of a polydispersity of rod lengths 
or a distribution of propulsion forces are not considered. The 
motion of rods are calculated with a discrete time step At ro d = 
0.001. Most of our rod simulations start from random initial 
states, where the rods are placed into the simulation box with 
random orientations and random positions without overlap. If 
not explicitly mentioned, the size of the simulation box is L x = 
L Y = 4004, which is much larger than the rod length. Periodic 
boundary conditions are employed. 

Our model differs from the model of Ref. l37ll by the type of 
repulsive interaction between the rods. In Ref. |[37l . rods in- 
teract by a "soft" volume exclusion, where the repulsion force 
is proportional to the square of overlapping area, while in our 
model the interaction is a short-range Lennard-Jones poten- 
tial between discretized beads. In the limit of a large overlap 
energy, the two models become equivalent. 



B. Flagella 

We consider a system of N/i flagella of length Ljt in a box 
of size L x x L y . Each flagellum consists of semi-flexible string 
of monomers of mass nifi, connected by springs (see Fig.[lJ)). 
The shape of the flagellum is determined by the elastic energy 



= X ^2 d R '-| " + Z ^3 { R '- +I " ^ocOR;} 2 + V . 



(9) 



Here, the first term is the harmonic potential generated by 
springs with spring constant k and rest length Iq. R, is the 
bond vector pointing from monomer i to monomer (z + 1). 
The second term of Eq. (|9]l is the bending energy the flag- 
ellum, with bending rigidity k. ^(Iqc) is an operator which 
rotates a two-dimensional vector clockwise by an angle Iqc. 
The local spontaneous curvature c varies with time t and posi- 
tion x along the flagellum to generate a propagating bending 
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velocities, 
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The detailed analysis of the beating pattern of nematodes 
and bull sperm 1 25 , 49 1 has shown that a single sine mode rep- 
resents the beating pattern to a good approximation. We use 
the wave number q = Itt/L/i, such that the phase difference 
between the first and the last monomer is 2n and one complete 
wavelength is present on the flagellum. The beating frequency 
/ is constant for each flagellum; it is chosen from a Gaussian 
distribution, centered at /o and with variance ct^/q . (p is the 
initial phase of the first monomer, which is chosen from a uni- 
form distribution in [0, 2n\. As t increases, a wave propagates 
along the flagellum from the first to the last monomer, push- 
ing the fluid backwards and propelling the flagellum forward. 
Although the spontaneous local curvature c is prescribed by 
Eq. ( 10 1, the flagellum is elastic and its configuration is af- 



fected by the viscosity of the medium and the flow field gen- 
erated by other flagella. The third term in Eq. (|9|l describes the 
interaction between flagella due to volume exclusion; here, we 
employ again the shifted and truncated Lennard-Jones poten- 
tial (Eq. (|6]l) between monomers of different flagella. 

Our model of a flagellum differs from the model of a sperm 
employed in Ref. [43 1 by the absence of a passive midpiece 
and a circular head. Also, in the sperm simulations ||43l , two 
sine waves were present on the tail, while a single sine wave 
is present on the flagellum. 

We use flagella of length Lfi = 50/q. The elastic moduli 
in Eq. ^ are the spring constant k — 1.25 x 10 4 k B T and the 
bending rigidity k = IQQksT Lji. The amplitude A = 5/Lfl of 



the spontaneous curvature in Eq. ( 10 1 induces a beating am- 
plitude of about 6.1 /o = 0.13 hf\. The strength e = 15ksT of 
the volume exclusion is large compared to the thermal energy. 
The simulations are initialized by placing Nfi flagella in the 
simulation box with random initial positions and orientations, 
without any overlap. The size of the simulation box is L x x L y , 
where L x = L Y = 4001q, eight times the length of a flagellum. 
Periodic boundary conditions are employed. 

Each simulation run of the flagella systems covers a total 
time interval of about 3300 beats. The first 800 beats are not 
taken into account in the calculation of averages, in order to 
allow the system to reach the stationary state. This time is 
longer than the largest relaxation time of about 650 beats ob- 
served in the system with a width <j>; = 0.1% of the frequency 
distribution. 



C. Multi-Particle-Collision Dynamics (MPC) 

MPC is a particle-based mesoscopic simulation technique 
used to describe the hydrodynamics and flow behavior of 
complex fluids. The fluid is modeled by N so i point particles 
of mass m so ij, which are characterized by their continuous 
space position r so y and velocity \ S oi,i- During every time step 
AtMPc, there are two distinct simulation steps, streaming and 
collision. In the streaming step, the fluid particles do not inter- 
act with each other and move ballistically according to their 



In the collision step, the particles are sorted into the cells of a 
square lattice of side length a according to their position, and 
interact with all other particles in same collision box through 
a multi-body collision. The collision step is defined by a ro- 
tation of all particle velocities in a box in a co-moving frame 
with its center of mass. Thus, the velocity of the 2-th particle 
in the j-th box after collision is 



\ soU (t + At MPC ) 
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where 



2 j Wlsol,iV solj 



(12) 



(13) 



is the center-of-mass velocity of y'-th box, and ¥lj(a) is a ro- 
tation matrix which rotates a vector by an angle +a, with the 
sign chosen at random. This implies that during the collisions 
particles exchange momentum, but the total momentum and 
kinetic energy are conserved within each collision box. In 
order to ensure Galilean invariance, a random shift of the col- 
lision grid has to be performed |50|. 

The total kinematic viscosity v is the sum of two contri- 
butions, the kinetic viscosity Vkm and the collision viscosity 
v co ii. In two dimension, approximate analytical expressions 
are||5l]|52, 
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where p is the average particle number in each box and h = 
AtMPc ^JksT lm so ia 2 is the rescaled mean free path. We use 
ksT = 1, m = 1, a = 1, AtMPC = 0.025, a = n/2, and 
p = 10. This implies, in particular, that the simulation time 
unit (ma 2 /ksT) 1 ^ 2 equals unity. With these parameters, the 
total kinematic viscosity of fluid is v = v co u + v^,„ « 3.02. 

During the MPC streaming step, the equations of motion 
of the flagella monomers are integrated using a velocity- 
Verlet algorithm, with a molecular-dynamics time step At/i = 
AtMPC /50 = 5xl0~ 4 . The bond length between the monomers 
is related to the collision cell size by lo = a/2. The flag- 
ella only interact with the fluid during the MPC collision step. 
This is done by sorting the flagella monomers together with 
the fluid particles into the collision cells and rotating their 
velocities relative to the center-of-mass velocity of each cell. 
Since energy is continuously injected into the system by the 
actively beating flagella, we employ a thermostat to keep the 
fluid temperature constant by rescaling all fluid-particle veloc- 
ities in a collision box relative to its center-of-mass velocity 
after each collision step. 

With the parameters given above, a single flagellum with 
fo = 1/120 swims forwards with the velocity v sing i e - 0.020 + 
0.001 in a MPC fluid. Thus, we estimate a Reynolds number 
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TABLE I: Definition of power-law exponents for the cluster size dis- 
tribution IT(;j), the average cluster size (/?), and the cluster lifetime 
Tiife, and their typical values for rods and flagella. 
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Tte = 2Af{V single I v » 0.04 for our flagellum model, where 
Afi = Q.12Lfi is the beating amplitude. The velocity of our 
flagella can be compared with the velocity of an infinitely long 
string beating in a two-dimensional fluid at Re = 0, which was 
calculated analytically by Taylor [24] to be 
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where A wave is the wave length and v wave - A wave f is the 
propagation velocity of the sinusoidal wave on the flagel- 
lum. Applying the parameters in our simulations, we obtain 
Vsingie = 0.0183, in excellent agreement with the simulation 
result. This demonstrates that the simulation model describes 
the limit of low -Reynolds-number hydrodynamics very well. 



FIG. 2: (Color) Snapshots of the SPR systems at different stationary 
states. Parameters are p r0 dL 2 md = 0.7744 and (a) 1/Pe = 0.02645 
(TIi); (b) 1/Pe = 0.00501 (TL.); (c) 1/Pe = 0.00041 (II3). Red dots 
mark the front ends of the rods, (d) Close-up of clusters of size n = 
3, 10 and 22 shows the partially blocked structure; chosen from a 
simulation with parameters p, dl^. od = 0.7744 and 1/Pe = 0.00095. 
For a movie see Ref. Ii53l . 



III. SWARMING BEHAVIOR OF SELF-PROPELLED 
RODS 



Cluster-Size Probability Density Functions and Stationary 
States 



After starting from a random initial state, the rods aggre- 
gate and form clusters. Large clusters can form by collisions 
of smaller ones, while at the same time they can break up due 
to collisions with other clusters or due to the noisy environ- 
ment. After a transient phase, the system reaches a stationary 
state, in which the formation rate of any cluster size equals 
its break-up rate. The degree of aggregation in the system 
depends on its parameters such as the Peclet number and the 
number density p ro j = N ro d/L x L y . 

We define a cluster as follows. We consider two rods to 
be in the same cluster if the angle between their orientation 
vectors is less than n/6 and the nearest distance is less than 
2lb, which is about two times the width of a rod. A cluster 
is defined as a set of rods that are neighbors either directly 
or through other rods at a given moment in time. Its size is 
simply the number of rods it contains. A freely gliding rod 
without any neighbor is considered as a cluster of size n = 1. 

We study systems at intermediate densities, where p m { is 
neither very low, such that there are hardly any collisions, nor 
high enough for a nematic phase to appear for rods in thermal 
equilibrium, i.e. densities lower than the critical density p c - 
3n/(2L? ,) of the isotropic-to-nematic phase transition l35ll . 



The statistical quantities, which will be analyze in Sees. Ill 
and|IV[ are listed in Tablejl] 



For a system with particles distributed at random in space, 
the probability of finding n particles in some area obeys a 
binomial distribution; in our SPR systems, the probability 
to find large particle numbers n is increased by aggregation 
and clustering. The stationary cluster-size probability den- 
sity function (PDF) II(n) results from the balance between the 
cluster formation and break-up rates. While the former de- 
pends on the collision rate of clusters, the latter depends also 
on the environmental noise and the cluster size. We distin- 
guish three different stationary states in our SPR systems by 
comparing the shapes of their corresponding PDFs. Snapshots 
are shown in Fig.|2j a movie can be found in Ref. [53 1. 

A disordered state, where rods are distributed in the whole 
space and oriented in different directions, is characterized by 
a PDF denoted as F^ in Fig. [3] In a snapshot (Fig. [2^), a weak 
aggregation tendency can be recognized in this case, where 
several small clusters of well polarized members glide in arbi- 
trary directions. FIi decreases as a power law for small cluster 
sizes, then decays exponentially for large n. The same kind of 
PDF has also been found in simulations of swarms of point- 
like SPP interacting via a phenomenological alignment mech- 
anism ifTTIfTBI . The range of the power-law-decay regime of 
FIj depends on the rod density and the environmental noise. 
Increasing density or decreasing noise shifts the exponential 
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FIG. 3: (Color online) Cluster-size distribution functions IT(«) for 
systems shown in the snapshots of Fig.|2ja), ITi (□, black), Fig.[2jb), 
n 2 (O, red), and Fig.gc), n 3 (A, blue). 



cut-off to larger n. 

The system with the second type of PDF, denoted II2 in 
Fig. [5] is more ordered, with an obvious tendency to form 
large clusters. A snapshot (Fig. ^p) shows several large and 
motile clusters moving in different directions. II2 also dis- 
plays a power-law decay at small cluster sizes, but shows an 
increased probability (compared to the power-law decay) of 
finding large clusters. Increasing the number density or de- 
creasing the noise shifts the prominent shoulder to larger clus- 
ter sizes. For very large aggregates, greater than the shoulder 
location, II2 decreases rapidly. 

The system with the third type of PDF, denoted II3 in Fig. [3] 
is characterized by a giant cluster, in which most rods are ori- 
ented radially towards the center (Fig. [2};). The giant cluster 
forms when several smaller motile clusters collide head-on in 
a short time interval, such that a nucleus with a blocked struc- 
ture emerges. This nucleus continues to grow until most of 
rods in the system are gathered in it. II3 has two parts, a peak 
at large n representing the giant clusters, and another peak at 
very small n corresponding to some freely swimming rods not 
collected by the giant cluster. The average rod density outside 
the giant clusters is very low. 

Both IT and n 2 display a power-law decay at small cluster 
sizes, 



(17) 



The exponent /3 is a function of the rod density p roc i and noise 
1/Pe; it increases with increasing p roi ] and decreases with in- 
creasing 1/Pe (Fig.Q. However, the dependence of /3 on p ro j 
or 1/Pe in the IT regime is much stronger than in the II2 
regime; in the latter case, /? approaches -2. 

By systematically varying the rod density p ro d and the envi- 
ronmental noise level, we can construct a phase diagram with 
regions characterized by different types of PDFs, see Fig. [5] 
Clearly, IT is found in the low-density and high-noise regime, 
II3 in the high-density and low-noise regime, and II2 is asso- 
ciated with the transition region between II 1 and II3. Note that 
all systems in Fig.[5]were started from disordered initial states. 
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FIG. 4: The exponent ft of the power-law part of IT and IT as a 
function of (a) the environmental noise 1 /Pe when prod^ B{l = 0.7744, 
and (b) the rod density p roi l} od when 1/Pe = 0.02645. 



Systems characterized by the probability density function II2 
bear some similarity with liquid systems supercooled below 
the freezing point. Note the system with 1/Pe — 0.00095 



and PmdL 2 rod — 0.11 '44 in Fig. [5jdisplays both II2 and II3 dis- 
tributions corresponding to simulations with different initial 
random states. Systems with the probability density function 
IJ3 show the characteristics of a glassy behavior, where the 
dense packing of rods arises from the random collisions, and 
remains frozen at later times. 

Our results are consistent with those of Ref. B71 . By com- 
paring short runs for systems with and without fluctuations, 
the transition from IT to II2 was found in Ref. l37l to shift to 
larger values of the aspect ratio L ro dlh and total area fraction 
of rods 7] = ProdLwdh- Fig- ^demonstrates that in our system 
the transition shifts with increasing 1/Pe to larger p ro dl? od , 
which is proportional to rjL ro( illb- 



B. Orientational Correlation Functions 

Although we distinguish three swarming states in our SPR 
systems, there are only two types of cluster structures. The 
motile clusters in the IT and II2 states consist of polarized 



7 




0.19 



0.38 



0.57 



0.76 



30 



■1.000 
— -0.7500 
\ -0.5000 
-0.2500 
I 
0.2500 
0.5000 
- 0.7500 
' 1.000 



FIG. 5: (Color online) Dynamical phase diagram of swarm behavior. 
Symbols indicate systems with PDF H[ (□, black), Yl 2 (O* rec D and 
II3 (A, blue). All systems were started from a random initial state. 
The dashed lines are guides to the eye. 

rods, as shown in Fig. [2^,b. In contrast, the giant clusters 
found in the II3 state consist of a large number of rods block- 
ing each other in their forward motion, as shown in Fig.[2j;. 

These two types of clusters can be distinguished by analyz- 
ing the orientational correlation function 

^"■ d^ ^*'''*-" 11 ' (18) 

Here u, is the unit vector denoting the orientation of rod i, 
ry(r, <p) is the vector pointing from the center of mass of rod 
i to rod j, and <p is the angle between u, and r,-,-. G(r) — > 1 
for r — > because two neighboring rods at close distance are 
always aligned. At large distance, G(r) — > 0. 

When the system is in a state characterized by II] or 1I2, 
G(r) is symmetric with respect to the direction = 0° with 
a maximum at r — (Fig. [6^). The slight elongation of G(r) 
in the directions <p — 0° and <p — 180° indicates that the clus- 
ters tend to slightly extend in the direction of the average rod 
orientation due to packaging. The width of G(r) is narrower 
in the front and wider in the back, because of their partially 
blocked structure (see Fig. |2}l) and because large clusters are 
more likely to collide with other clusters head-on. If an head- 
to-head collision does not result in the formation a larger clus- 
ter or a blocked structure, the front tips are sharpened due to 
the "attrition" of the two clusters. 

If the system is in the state with a giant cluster, G(r) shows 
a very different behavior, see Fig. [6j?. G(r) still has a positive 
maximum near r — 0, which represents a high local orien- 
tational order. However, a region with negative correlations, 
G < 0, develops, with a minimum at some (/ ,<p'). Because all 
rods point preferentially towards the center of cluster, the pro- 
pelling forces of the rod nearly compensate each other. There- 
fore, the locomotion speed of a giant cluster is much smaller 
than the gliding speed of a single rod. Moreover, the pro- 
pelling forces generate a net torque due to the deviation of the 
rod orientations from pointing exactly towards the center of 
mass, which implies a rotational motion of the giant cluster. 




240 ' 300 

270 



FIG. 6: (Color) The orientational correlation function G(r) as a func- 
tion of the relative position r (a) in a system with IT (prodL^ = 
0.7744 and 1/Pe = 0.00501), where there are motile clusters; (b) in 
a system with II3 (p rod l} rod = 0.7744 and 1/Pe = 0.00041) character- 
ized by the presence of a giant cluster. 



<f>' is related to this rotation. For 0° < <f>' < 90°, the cluster 
rotates counterclockwise; for -90° <(/>'< 0°, it rotates clock- 
wise; for <f>' = 0°, there is no net torque and the giant cluster 
does not rotate. 



C. Average Cluster Size 

The average cluster size («) of the system is 

<n>=^«n(n), (19) 

n 

where Yl(n) is the normalized cluster-size distribution func- 
tion. («) increases with increasing p roc /, as shown in Fig. [7^; 
in the low-density limit, («) approaches unity. («) decreases 
with increasing noise level, 1/Pe, as shown in Fig.[7j3. In the 
TI2 regime, the system exists in two metastable states, depend- 
ing on the initial conditions. With random initial conditions, 
a "supercooled" state emerges, which transforms into the n3 
state once a giant-cluster nucleus has formed. This can be seen 
in Fig. [T]} for 1 /Pe = 0.00095, where two data points show 
simulation results with different random number for random 
initial states. With a giant cluster as initial state, the system 
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FIG. 7: (Color online) The average cluster size (n) as a function of 
(a) the rod density ProdL 2 rod with 1/Pe = 0.02645 and (b) the envi- 
ronmental noise 1 /Pe with p roi Lr rod = 0.7744. The solid (red) line 
is a fit of the power-law part with exponent C, = 0.37. The dashed 
(black) lines are the boundaries separating different regions in the 
phase diagram (Fig.|5j. The open symbols represent systems with 
IT (A, blue), IT2 (0> rec 0, and IT (□, black), starting from random 
initial states. The solid symbols represent the systems with n 2 (a, 
blue) and IT (■, black), starting from a state with a giant cluster. 



stays in the II3 state unless the noise is large enough to de- 
stroy the giant cluster; this occurs in Fig.|7j3 for 1 /Pe = 0.04. 
Interestingly, (n) shows a power-law decay 



<«> ~ Pe f 



(20) 



in the 111 and II2 region when the system starts from a disor- 
dered state, with exponent £ 0.37. 



D. Cluster Lifetime 



FIG. 8: (Color online) Average cluster lifetime 7jf/ e (n) for systems 
with the same rod number density p roi l2 oi = 0.7744 but with a dif- 
ferent noise level, as indicated. The dashed lines are power laws l |21| l 
with an exponent S = 0.2. 



single-rods "clusters" cannot disintegrate. For n > 2, Tuf e (n) 
decreases smoothly with increasing cluster size. The data for 
mid-size clusters (2 < n < 30) show an effective power-law 
dependence, 



Tufe(n) 



(21) 



with an exponent 6 0.2. Because the environmental noise 
determines the break-up rate of clusters, Tuje increases with 
decreasing 1/Pe. We only show the lifetime of motile clusters 
in systems characterized by lli and II2. The giant clusters 
found in the state characterized by II3 can persist for a very 
long time until a sufficiently large fluctuation occurs. 

To understand the dependence of the cluster lifetime on n, 
we can assume that only single rods are lost at the cluster sur- 
face ll37ll . In this case, the probability to loose a rod per unit 
time is proportional to the perimeter length, which scales as 
n 1 ^ 2 (for compact clusters of approximately circular shape). 



Therefore, this simple argument implies a scaling law (21 
with exponent 6 = 0.5. The growth of clusters is more com- 
plex, since it can occur by collision with all types of other 
clusters; however, the collision cross-section should again be 
proportional to n 1 ' 2 . The value of 5 — 0.5 is considerable 
smaller (corresponding to shorter lifetimes for larger clusters) 
than observed in our simulations. This indicates that there 
must be another mechanism of cluster decay. Indeed, the typ- 
ical cluster configurations of Fig. |2jl indicate that only at few 
places along the perimeter, rods may have the possibility to 
leave the cluster. 



We define the lifetime of a cluster as the length of the time 
during which its members do not change. The lifetime of a 
cluster is analyzed with a time interval At = 100; thus, cluster 
lifetimes less than At cannot be resolved. The average cluster 
lifetime Tu/ e is a function of cluster size n. 

As shown in Fig. [8] the lifetimes of the clusters of size n — 1 
are always much longer than of other cluster sizes, because 



E. Finite-Size Effects 

In our simulations, the finite simulation-box size implies a 
finite number of particles. A cluster can never grow larger 
than the total number of rods in the system. Consequently, 
all quantities related to the cluster size, such as the cluster 
size distribution II and the stationary average cluster size («) 
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FIG. 9: (Color online) Effect of finite system size on the probability 
density distribution function Ht(n) for p md L 2 md = 0.7744, 1/Pe = 
0.04009, and different simulation box sizes, as indicated. The inset 
shows the exponent ft of the power-law part of H\(n) as a function of 
the size L x of the simulation box. 



display finite-size effects. Similarly, density fluctuations at the 
scale of the simulation-box size are suppressed. 

For the probability density function El(n), the absence of 
cluster larger than N roc i does not only introduce a cut-off at 
large cluster size, but also affects the exponent ji of the power- 
law part, as shown in Fig. [9] For systems with Fli, the data for 
small box sizes (L x = L y = 4.5L ro d and 13L ro d) still obey a 
power-law decay at small n, without an obvious change of the 
exponent, as shown in the inset of Fig.[9} but they deviate from 
the power law when n approaches N ro d- When the simulation 
box is large (L x - L y - 36Zw and 51L, j), the PDFs almost 
coincide, and their exponential cut-offs are observed at the 
same value of n; also, /3 approaches an asymptotic value when 
L x increases. Therefore, we conclude that our results for the 
larger systems represent the thermodynamic limit. Similarly, 
the power-law part of II2 extends with increasing box size, and 
the location of the prominent shoulder shifts to larger cluster 
size. The finite-size effects are significantly stronger for sys- 
tems in the II3 region of the phase diagram. When the system 
is too small, the total rod number is not sufficient to trigger 
the formation of a blocked structure. The system then stays 
in a IT2 state. This supports the claim that the state with IT2 is 
a "supercooled" state. We believe that the absence of the IT3 
state in Ref. 11371 is due to finite-size effect; a system of only 
100 rods is too small to form a blocked structure. 

The dependence of the average cluster size («) on the lin- 



ear system size L x is shown in Fig. 10 For systems with EE 
and 1I2, («) increases with L x and eventually reaches a plateau 
value. For the system with EI3, («) strongly diverges when L x 
increases. Thus, («) can be considered as an intensive quan- 
tity in the first two states, and as an extensive quantity in the 
third state. 

Suppose the probability density function EI(n) obeys a 
power law for all cluster sizes, 



n(«) 
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FIG. 10: (Color online) The average cluster size (n) as a function of 
the size L x of the simulation box for the three clustering states. The 
number density is p ro dL 2 rod = 0.7744 in all systems. 



where p < — 1 and N = p ro dL x L y is the total number of rods 
in the system. It is easy to verify that for N » 1, where sums 

rN 

over n can be well approximated by integrals, J { Tl(n)dn = 1, 
so that Il(n) is properly normalized. The sharp drop due to the 
limited box size is neglected. In this case, the average cluster 
size of the system is obtained to be 



(n) 



' -(1 + /3)N 2+l3 /(2 +j3), -2<p<-\ 

(1 -ArVln/y, jS = -2 
. (l +y e)/(2+j6), yS<-2 



(23) 



For -2 < /3 < -1, the average cluster size strongly depends 
on the total number of rods, whereas for /3 < -2, («} is 
independent of A^. For large negative f3, (n) approaches unity, 
which means that all rods are gliding freely. 

In our simulations, the effective exponents in the and IT2 
regimes are -6 < < -2.5 and -2.5 < f3 < -2.0, respectively, 
see Fig. [4] Thus, Eq. ( 23 1 implies that finite-size effects are 
weak in the Fli regime, and are pronounced in the IT2 regime, 



in agreement with the simulation results of Fig. 10 (Eq. (23 1 



does not apply to the II3 state since the assumption of a power- 



law dependence (22 1 does not hold.) 



IV. 



SWARMING BEHAVIOR OF FLAGELLA IN A MPC 
FLUID 



Multi-flagellum systems show a similar swarming behav- 
ior, consisting of aggregation and clustering, as observed in 



(22) 



Sec. Ill for self-propelled rods (see Fig. 11 and movie [53]). 
Synchronization of the flagellar beat, and attraction and align- 
ment of flagella do not only arise from volume exclusion, as 
in the SPR systems, but are also triggered by the hydrody- 
namic interactions between the sinusoidally undulating bod- 
ies l43ll44l . At the same time, hydrodynamic interactions be- 
tween flagella act as a source of environmental noise, which 
causes the flagella trajectories to fluctuate strongly. 
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FIG. 1 1 : (Color online) Snapshot of a multi-flagellum system in a 
MPC fluid. The parameters are P//Ly 7 = 1.5625 and crp = 0.1%. 
The black box shows the boundary of the simulation box. Periodic 
boundary conditions are employed. For a movie see Ref. 1531 . 




FIG. 12: (Color online) Synchronization and attraction of two fiag- 
ella. The flagella have the same beating frequency / = 1/120, and a 
phase difference Atp = 0.5n. The snapshots are taken at the times (a) 
tf = 0.167, (b) tf = All and (c) tf = 22.2. 



A. Hydrodynamic Synchronization, Attraction, and 
Aggregation 

The synchronization and attraction of two flagella is shown 
in Fig. 



12 



Synchronization is achieved within about four 
beats, while the formation of a tight pair from an initial dis- 
tance of about one-third of the flagellar length takes about 
20 beats. The flow field of a flagellum is shown in Fig. 13 



The flow field at a certain time in the beating cycle (Fig. 13 1) 



shows that formation of two vortices, which propagate from 
the front to the rear end as the flagellum moves forward. 

The hydrodynamic interaction of swimmers depends on the 
type of self-propulsion. The average flow field of flagellum, 
integrated over the whole beating cycle, demonstrates that the 
flagellum, which might be expected to be a "neutral" swim- 
mer (i.e., neither a pusher nor a puller) is indeed a very weak 
pusher — where the dominant propulsion is located closer to 
the rear end — because the line connecting the centers of the 
two vortices intersects the average flagellum shape behind its 



mid-point (Fig. 13 ~>). This generates a in-flow from both sides 



of the flagellum near the front end, which is responsible for 
hydrodynamic attraction fiOlPffl . 

In multi-flagellum systems, large clusters can form by col- 
lisions of smaller clusters, supported by the hydrodynamic at- 
traction between neighboring flagella; large clusters can disin- 
tegrate into smaller components due to the diversity of flagel- 
lar frequencies, or the hydrodynamic flow fields of other clus- 
ters. With hydrodynamic interactions, large clusters of flag- 
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FIG. 13: (Color online) Flow field of a single flagellum, (a) at a fixed 
time in the beating cycle, and (b) averaged over the whole beating 
cycle. A snapshot of the flagellum and the average flagellum shape 
are superimposed in (a) and (b), respectively. The scale bar indicates 
the magnitude of the flow velocities. 



ella are usually strongly extended in their direction of motion, 
as shown in Fig. [TT]and movie |[53ll . The flagella inside the 
cluster are well synchronized. This structure is reminiscent of 
the "sperm-train" structure observed in rodent-sperm experi- 
ments [27, 28 1 . The elongated clusters can extend to distances 
as large as the side length of the simulation box, which in- 
duces strong finite-size effects. 



Similar to the definition of a rod cluster in Sec. Ill a flag- 
ellum cluster is defined as a set of flagella that are connected 
or neighbors either directly or through other agents at a given 
moment in time. Its size is the number n of flagella it contains. 
A freely-swimming single flagellum is considered as a cluster 
of size n = 1. 



B. Cluster-Size Distributions 



Both probability density functions II! and FI2 are observed 



in our multi-flagellum systems, as shown in Fig. 14 The vari- 



ance <Tfi of the distribution of beat frequencies is used as a 
measure of the noise level. At low pp or high cry/, we find FTi; 
at high pfi or low cry?, we observe II2. In contrast to II2 for 
SPR systems in Sec. Ill A II2 for flagella systems displays a 
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FIG. 14: (Color online) The two different types of cluster-size prob- 
ability density functions Ht(n) with cry; = 2% (□, black), and n 2 (rc) 
with 07; = 0.1% (O, red), observed in multi-flagellum systems with 
density PfiL 2 fl = 1.5625. The dashed lines are fits to the power-law 
parts of each PDF. The inset shows the exponent ft as a function of 
the variance cry; of the frequency distribution. Symbols indicate the 
systems characterized by IT (0> red) and IT (□, black). The dashed 
line is a guide to the eye. 




FIG. 15: (Color online) Stationary average cluster size (n) versus 
the variance (Tfi of the frequency distribution, for flagellar density 
P/iL^j = 1.5625. Symbols indicate systems with IT (□, black) and 
II2 (0> re d). The error bars are the standard deviation of the fluctua- 
tions. The dashed line indicates a power-law decay with an exponent 
( = 0.26. 



deviation from the power-law behavior for very small cluster 
sizes, n — 1 and n — 2. We believe that this is due to the 
hydrodynamic synchronization and attraction of neighboring 
fiagella. For flagella, we have never observed a giant clus- 
ter with a blocked structure, in contrast to the SPR system of 
Fig. 12};. 

Although the distribution of beating frequencies is an in- 
ternal property of the swimmers, the influence of cry; on the 
exponent ft of Eq. ( fT7| ) is similar to the influence of the envi- 
ronmental noise in our previous SPR simulations, as shown in 
the inset of Fig. [14] /3 is nearly constant for cry; < 3%, then 
decreases smoothly with increasing cry/. 

The average cluster size (n) in the stationary state is a func- 



FIG. 16: (Color online) The lifetime Tnf e v sing i e ILfi of flagella clus- 
ters versus their size n. The flagella number density is P//Ly; = 
1 .5625. The dashed lines indicate a power law with exponent 6 = 0.5. 



tion of (Tfi, as shown in Fig. 15 Increasing cry; results in an in- 
crease of the overall break-up rate; hence (ri) decreases. In the 
large cry; limit, (n) approaches unity, corresponding to a disor- 
dered state with randomly distributed flagella. The power-law 
decay 



<«> 



'fi 



(24) 



of the average cluster size with exponent ( =t 0.26 emphasizes 
the universality of the swarming behavior of rSPP systems in 
two dimensions. The power-law scaling of (n) as a function of 
(jfi implies a divergence when cry; — > 0. We believe that the 
small deviation from the power-law behavior for cry; = 0.1% 
in Fig. 15 as well as the deviation of f3 from the plateau value 
for (Tfi = 0.1% in Fig. 14 are due to finite-size effects. 



C. Cluster Lifetimes 



The average cluster lifetime Tnf e (n) decreases as an effec 



tive power-law function of cluster size «, see Eq. (21 1, with an 



exponent 6 0.5, as shown in Fig. 16 The value of 6 is very 
close to the prediction based on the assumption of a mech- 
anism of particle accumulation and shedding proportional to 
the cluster perimeter, as presented in Sec. HID This good 
agreement provides further evidence for the different mecha- 
nisms of cluster stabilization for rods and flagella, which are 
a (partially) blocked motion and a hydrodynamic attraction, 
respectively. 

Note that the system size of the flagella simulations is not 
as large as for the SPR systems. Thus, the effective power law 
can only be observed over a smaller range of cluster sizes. In 
SPR simulations, single rods (n — 1) always have a much 
longer lifetime compared to expectation from the effective 
power law, see Fig. [8] In contrast, for flagella with full hy- 
drodynamic interactions, Tnf e (l) is much closer to the power- 
law extrapolation, and can even be lower than the power-law 
prediction (e.g. for cryy = 0.5% in Fig. 16 1. 
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FIG. 17: (Color online) Flow field of a single sperm, beating with 
two sine waves on its tail, averaged over the whole beating cycle. 
The average sperm shape is superimposed. The scale bar indicates 
the magnitude of the flow velocities. 



D. Comparison of Sperm and Flagella 



As explained in Sec. II B 



our model of a fiagellum differs 
from the model of a sperm employed in Ref. [43 1 by the ab- 
sence of a passive midpiece and a circular head. Also, in the 
sperm simulations 11431 . two sine waves were present on the 
tail, while a single sine wave is present on the fiagellum. 

How similar or different is the collective behavior of sperm 
and flagella? There are three different aspects to this ques- 
tion. Synchronisation depends mainly on the interaction of 
the time-dependent oscillatory flow field of two neighboring 
flagella G4ll44l . and is therefore very similar, as can be seen 



from the results presented in Sec. IV A and those of Ref. [43|. 
On the other hand, the hydrodynamic attraction of sperm and 
flagella is quite different. A sperm cell, consisting of a fiag- 
ellum and a large head, is clearly a pusher, as demonstrated 
by the average flow field of a sperm in Fig. 17 The fiagellum 



pushes the fluid backward in both cases, but the bulky head of 
the sperm drags the fluid forward much more strongly, which 
generates the characteristic sidewise inflow of fluid towards 
the midpiece region |40, 4T)|45]. In contrast, flagella are very 



weak pushers, as demonstrated in Fig. 13 *> above. Therefore, 
sperm have stronger hydrodynamic attraction than flagella. 

Finally, the swarming behavior in both flagella and sperm 
system is characterized by cluster-size distributions and the 
dependence of the average cluster size on the width cr/i of the 
distribution of beat frequencies. While the cluster-size dis- 
tribution of flagella follows a power-law decay over a wide 
range, it was not possible to clearly identify a power-law be- 
havior for sperm in Ref. [43 1 due to the relatively small sys- 
tems of 25 and 50 sperm. The average cluster size is found 
to depend on <x// as (n) ~ o-Jj, with £ = 0.20 for sperm [43 1 
and £ = 0.26 for flagella. Larger systems have to be inves- 
tigated to see whether the exponents £ for sperm and flagella 
are different or not. In any case, the stronger hydrodynamic 
attraction of sperm, which favors larger cluster sizes, is par- 
tially offset by the bulky head of sperm, which implies that 
the sperm clusters in Ref. [43 1 are much more loosely packed 



than the flagella clusters studied here. 



V. SUMMARY AND CONCLUSIONS 

We have simulated systems of rigid rods propelled by a 
constant force along their long axis, and systems of flagella 
propelled by a sinusoidal beating motion, in two dimensions. 
In both systems, we observe cluster formation and break-up, 
controlled by the particle density and the internal or external 
noise. In our simulations, the particle density is always much 
lower than the critical density of a nematic phase in thermal 
equilibrium. 

Without any attractive potential, self-propelled rods (SPR) 
exhibit an aggregation behavior triggered only by volume ex- 
clusion. Three characteristic types of cluster-size probabil- 
ity density functions IT(n) appear in different regions of a dy- 
namic phase diagram of stationary states. At high noise and 
low density, the system is characterized by LTi, which shows 
a power-law distribution over a range of cluster sizes, with an 
exponential cutoff at large cluster sizes. At low noise and high 
density, the system is in a state characterized by IT3, which has 
a peak at sizes near the total number of particles in the system, 
representing a giant cluster. Systems in an intermediate region 
of noise and density are characterized by IL, which is a tran- 
sition state between LIi and IT3. It has a bimodal shape, with a 
power-law decay at small cluster sizes and a shoulder at larger 
sizes. Clusters in IF and FI2 systems retain a high motility, 
whereas the giant clusters found in the third state is almost 
immobile due to its blocked configuration. The average clus- 
ter size at equilibrium, directly related to cluster-size distri- 
bution II, displays a power-law dependence with decreasing 
noise amplitude before the system reaches the II3 state. 

Sinusoidally beating flagella were simulated in a low- 
Reynolds -number fluid with full hydrodynamics as an exam- 
ple of self-propelled rod-like particles with explicit propulsion 
mechanism. Flagella synchronize their beats and attract each 
other through the hydrodynamic interactions. Despite the dif- 
ferent propulsion mechanisms, the basic swarm behavior of 
aggregation and clustering observed for swimming flagella is 
remarkably similar to the behavior seen in SPR systems. We 
observe both IF and IT2 cluster-size probability density func- 
tions by varying the width cry; of the flagellar beat-frequency 
distribution, which acts as a source of internal noise in the 
system. The average cluster size also display a power-law de- 
pendence on crfi, as for SPR systems. 

Despite these similarities in the clustering behavior, the 
two systems show some important differences. They can be 
traced back to the hydrodynamic attraction between beating 
flagella, which is absent in our simulations of self-propelled 
rods. First, the configurations of the flagella clusters consist 
of tightly stacked flagella with synchronized shapes, and ex- 
tend in their moving directions. Those elongated clusters are 
reminiscent of the huge, mobile "sperm trains" observed in 
rodent-sperm experiments [27|. Clusters in the SPR systems 
are more compact, and have a wedge-like structure, which 
arises from the partially blocked rod motion responsible for 
the cluster aggregation, as well as from collisions with other 
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clusters. Second, the II3 state of a completely blocked struc- 
ture, which is observed for SPR at high density and low noise, 
does not seem to exist in flagellar systems. Third, the cluster 
lifetimes decay with different effective power laws, 5 = 0.2 
for SPR and 6 = 0.5 for fiagella. Finally, hydrodynamic inter- 
actions between different fiagella clusters act as an additional 
source of noise and contribute to increase the break-up rate. 

The existence of the giant, immobile cluster should depend 
sensitively on the aspect ratio and the type and range of the in- 
teractions between self-propelled rods, where longer rods and 
shorter-range interaction favors the giant-cluster formation. 
This conclusion follows from the result of Ref. [ 37] for rods of 
aspect ratio L ro d/lb < 12 that IT-ITi boundary shifts to higher 
density with decreasing rod length, and our result of Fig.[7jthat 
the II2 state corresponds to "supercooled" liquid state which 
transforms into the TI3 state once a giant-cluster nucleus has 
formed. Blocked clusters were not seen in Ref. 071 for rod 
lengths L ro dllb < 12 due to the relatively small system size 
with N roc i = 100. However, blocked states were observed in 
Ref. 1T361 for a much larger rod length, L ro d/lb = 40, already 
for a system of only about 50 rods at density p rol jL? od = 2. 

Our simulations have been restricted to the isotropic phase 



of rods in thermal equilibrium. It will be interesting to see 
in the future whether immobile, blocked states can also exist 
(or even dominate) in the nematic regime, or whether they are 
suppressed by the preferred rod orientation. 

In the light of our results, we conclude that different sys- 
tems of rod-like self-propelled particles display a universal 
swarming behavior, but also specific properties related to their 
propulsion mechanisms and the presence or absence of hydro- 
dynamic interactions. 
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